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Abstract The generation of pseudo-random discrete 
probability distributions is of paramount importance 
for a wide range of stochastic simulations spanning from 
Monte Carlo methods to the random sampling of quan¬ 
tum states for investigations in quantum information 
science. In spite of its significance, a thorough expo¬ 
sition of such a procedure is lacking in the literature. 
In this article we present relevant details concerning 
the numerical implementation and applicability of what 
we call the iid, normalization, and trigonometric meth¬ 
ods for generating an unbiased probability vector p = 
(pi, • • • ,Pd)- An immediate application of these results 
regarding the generation of pseudo-random pure quan¬ 
tum states is also described. 

Keywords Discrete Probability Distributions • 
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1 Introduction 

Roughly speaking, randomness is the fact that, even 
using all the information that we have about a physical 
system, in some situations it is impossible, or unfeasi¬ 
ble, for us to predict exactly what will be the future 
state of that system. Randomness is a facet of nature 
that is ubiquitous and very influential in our and other 
societies As a consequence, it is also an essen¬ 

tial aspect of our science and technology. The related 
research theme, that was motivated initially mainly by 
gambling and led eventually to probability theory [4l 
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[^, is nowadays a crucial part of many different fields 
of study such as computational simulations, informa¬ 
tion theory, cryptography, statistical estimation, sys¬ 
tem identification, and many others [HllTlIMI^ ITOlfTniT^ . 

One rapid-growing area of research for which ran¬ 
domness is a key concept is the maturing field of quan¬ 
tum information science (QIS). Our main aim in this 
interdisciplinary field is understanding how quantum 
systems can be harnessed in order to use all Nature’s 
potentialities for information storage, processing, trans¬ 
mission, and protection [nmna. 

Quantum mechanics [Bin] is one of the present 
fundamental theories of Nature. The essential mathe¬ 
matical object in this theory is the density operator (or 
density matrix) p. It embodies all our knowledge about 
the preparation of the system, i.e., about its state. From 
the mathematical point of view, a density matrix is sim¬ 
ply a positive semi-definite matrix (notation: p > 0) 
with trace equal to one (Tr(p) = 1). Such kind of ma¬ 
trix can be written as p = which is known as 

the spectral decomposition of p. In the last equation 11 j 
is the projector {Ujllk = Sjkllj and = Id, where 

Id is the d—dimensional identity matrix) on the vector 
subspace corresponding to the eigenvalue Vj of p. From 
the positivity of p follows that, besides it being Hermi- 
tian and hence having real eigenvalues, its eigenvalues 
are also non-negative, rj > 0. Once the trace function 
is base independent, the eigenvalues of p must sum up 
to one, Tr(p) = = 1- Thus we see that the set 

{rj} possesses all the features that define a probability 
distribution (see e.g. Ref. HI)- 

The generation of pseudo-random quantum states is 
an essential tool for inquires in QIS (see e.g. Refs. [B 

and involves two 

parts. The first one is the generation of pseudo-random 
sets of projectors {Uj}, that can be cast in terms of the 
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creation of pseudo-random unitary matrices. There are 
several methods for accomplishing this last task 
mm, whose details shall not be discussed here. Here 
we will address the second part, which is the generation 
of pseudo-random discrete probability distributions [85l 
[3132], dubbed here as pseudo-random probability vec¬ 
tors (pRPV). 

In this article we go into the details of three methods 
for generating numerically pRPV. We present the prob¬ 
lem details in Sec. m The Sec. |3I is devoted to present 
a simple method, the iid method, and to show that it 
is not useful for the task regarded in this article. In 
Sec. [4] the standard normalization method is discussed. 
The bias of the pRPV appearing in its naive-direct im¬ 
plementation is highlighted. A simple solution to this 
problem via random shuffling of the pRPV components 
is then presented. In Sec.jSjwe consider the trigonomet¬ 
ric method. After discussing some issues regarding its 
biasing and numerical implementation, we study and 
compare the probability distribution generated and the 
computer time required by the last two methods when 
the dimension of the pRPV is varied. The conclusions 
and prospects are presented in Sec. |6| 

2 The problem 

By definition, a discrete probability distribution [3] is 
a set of non-negative real numbers, 

Pj > 0, (1) 

that sum up to one, 

EUp^ = 1 - ( 2 ) 

In this article we will utilize the numbers pj as the 
components of a probability vector 

p= {pi,--- ,Pd)- (3) 

Despite the nonexistence of consensus regarding the 
meaning of probabilities [4], here we can consider pj 
simply as the relative frequency with which a particu¬ 
lar value Xj of a physical observable modeled by a ran¬ 
dom variable X is obtained in measurements of that 
observable under appropriate conditions. 

We would like to generate numerically a pseudo¬ 
random probability vector p whose components {pj}'j^i 
form a probability distribution, i.e., respect Eqs. o 
and (ED- In addition we would like the pRPV to be 
unbiased, i.e., the components of p must have similar 
probability distributions. A necessary condition for ful¬ 
filling this last requisite is that the average value of pj 
(notation: (pj)) becomes closer to 1/d as the number of 
pRPV generated becomes large. 


At the outset we will need a pseudo-random number 
generator (pRNG). In this article we use the Mersenne 
Twister pRNG [38], that yields pseudo-random num¬ 
bers (pRN) with uniform distribution in the interval 
[0,1]. We observe however that the results reported in 
article can be applied straightforwardly also when deal¬ 
ing with true random numbers 1311110]. 


3 The iid method 

A simple way to generate an unbiased pseudo-random 
probability vector p = (pi, • ■ • ,Pd) is as follows. If we 
create d independent pseudo-random numbers Xj with 
identical probability distributions in the interval [0,1] 
(so the name of the method) and set 


we will obtain a well defined discrete probability distri¬ 
bution, i.e., Pj > 0 and E‘j=iPj — 1- Besides, as p is 
unbiased, the mean value of pj approaches 1/d as the 
number of samples grows. 

Nevertheless, we should note that the sum Ek=i 
shall be typically greater than one. This in turn will 
lead to the impossibility for the occurrence of pRPVs 
with one of its components equal (or even close) to one. 
As can be seem in Fig. El this problem becomes more 
and more important as d increases. Therefore, this kind 
of drawback totally precludes the applicability of the iid 
method for the task regarded in this article. 
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Fig. 1 (color online) Probability distribution for the 
first component of the unbiased probability vector p = 
(pi,--- ,pd) for one million random samples generated us¬ 
ing the iid method. In the inset is shown the probability 
distribution for the components of the probability vector 
P = (Pl,P2,P3,P4). 
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4 The normalization method 

Lets us begin our discussion of this method by consid¬ 
ering a probability vector with dimension d = 2, i.e., 
p = {pi,p 2 )- If the pRNG is used to obtain pi G [0,1] 
and we impose the normalization to get p 2 = 1 — pi, we 
are guaranteed to generate an uniform distribution for 
Pj G [0,1] for both j = 1 and j = 2. 

If d = 3 then p = {pi,P 2 ,P 3 ) and the pRNG is 
used again (two times) to obtain pi G [0,1] and p 2 G 
[0,1 — pi\. Note that the interval for p 2 was changed 
because of the normalization of the probability distri¬ 
bution, which is also used to write Pa = I — (pi + P 2 )- 
As Pi is equiprobable in [0,1], for a large number of 
samples of the pRPV, its mean value will be 1/2. This 
shall restrict the values of the other components of p, 
shifting the “center” of their probability distributions to 
1/4, biasing thus the pRPV. Of course, if one increases 
the dimension of the pRPV, the same effect continues 
to be observed, as is illustrated in the table below for 
10® pRPV generated for each value of d. 


d 

(Pi) 

(P 2 ) 

(Pa) 

(P4) 

(Ps) 

2 

0.5001 

0.4999 




3 

0.5000 

0.2499 

0.2501 



4 

0.4999 

0.2503 

0.1248 

0.1250 


5 

0.4998 

0.2501 

0.1252 

0.0625 

0.0624 


The probability distributions for the four compo¬ 
nents of the probability vector p = (pi,P 2 ,P 3 ,P 4 ) are 
shown in Fig. [2] For the sake of illustration, the spaces 
for the pRPV are sketched geometrically in Fig. |3] for 
the cases d = 2 and d = 3. We observe that the proce¬ 
dure for generating p as explained above is the moti- 



Fig. 2 (color online) Probability distribution for the compo¬ 
nents of the biased probability vector p = (pi,P 2 ,P 3 ,P 4 ) for 
one million random samples of it. In the inset is shown the 
probability distribution for the components of the unbiased 
probability vector q = (gi, 92 , 93 , 94 )■ 




Fig. 3 (color online) In blue is shown the set of possible 
values for the components of the probability vector in the 
cases of d = 2 (figure on the top) and for d = 3 (figure on 
the bottom). For higher dimensions the probability space is 
a hyperplane. 


vation for the name of the method, the normalization 
method. 

A simple solution for the biasing problem just dis¬ 
cussed is shuffling the components of the pRPV in each 
run of the numerical experiment. This can be done, 
for example, by generating a random permutation of 
{1, 2, • • • , d — 1, d}, let us call it {ki,k 2 , • • • , kd-i,kd}, 
and dehning a new pRPV as 

q = ,qd-i,qd) 

■■= {Pki,Pk2,--- ,Pk^-l,Pk^)- (5) 

In the table below is presented the mean value of the 
components of q for 10® pRPV generated. 


d 

(91) 

( 92 ) 

( 93 ) 

( 94 ) 

( 95 ) 

2 

0.5005 

0.4995 




3 

0.3332 

0.3331 

0.3337 



4 

0.2494 

0.2503 

0.2499 

0.2504 


5 

0.2001 

0.1997 

0.2006 

0.1999 

0.1997 


In the inset of Fig. [2] is shown an example with the 
resulting probability distributions for the four compo¬ 
nents of q = ( 91 , 92 , 93 , 94 )- 

From the discussion above we see that in addition 
to the d — 1 pRN needed for the biased pRPV, we have 
to generate more d — 1 pRN for the shuffling used in 
order to obtain an unbiased pRPV (because 3 — 
d{d + l)/2), resulting in a total of 2{d — 1) pRN per 
pRPV. 
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5 The trigonometric method 

As the name indicates, this method uses a trigonomet¬ 
ric parametrization for the components of the 

probability vector p = (pi, • • • ,pd)- 

d-l 

Pj := sin^ 0j-i cos^ 9k, (6) 

k=j 

with 9q = 7r/2 (so the name trigonometric method). 
More explicitly, 

cy r\ r\ r\ cy 

Pi = sin 00 cos 9 1 cos 02 cos 03 • • • cos 9d-i 
P2 = sin^ 01 cos^ 02 cos^ 03 • • • cos^ 9d-i 
P3 = sin^ 02 COS^ 03 COS^ 04 • • • COS^ 9d-l 

Pd-1 = sin^ 9d-2 COS^ 9d-i 


more terms to have values closer to zero. Thus, there is 
the issue of a biased pRPV again. A possible solution 
for this problem is, once more, shuffling. In the next two 
tables are shown the average value of the components 
of 10® pRPV generated via the trigonometric method 
before, (pj), and after, (qj), shuffling. 


d 

(pi) 

(P2) 

(ps) 

(Pi) 

(Ps) 

2 

0.4999 

0.5001 




3 

0.2502 

0.2498 

0.4999 



4 

0.1250 

0.1251 

0.2497 

0.5003 


5 

0.0625 

0.0624 

0.1250 

0.2496 

0.5006 


d 

( 01 ) 

( 02 ) 

( 03 ) 

( 04 ) 

( 05 ) 

2 

0.5006 

0.4994 




3 

0.3331 

0.3334 

0.3337 



4 

0.2502 

0.2497 

0.2501 

0.2450 


5 

0.2005 

0.1998 

0.1999 

0.2001 

0.1997 


Pd = sin^0d_i. (7) 

A simple application of the equality cos^ 9j + sin^ 9j = 1 
to this last equation shows that pj > 0 and J2j=i Pj — 
1. Therefore this parametrization, which utilizes d — 1 
angles 9j , leads to a well defined probability distribution 
{Pj}j=i- 

Let us regard the numerical generation of an unbi¬ 
ased pseudo-random probability vector by starting with 
the parametrization in Eq. ([S]). Of course, this task is 
accomplished if the components of the pRPV are cre¬ 
ated with uniform probability distributions. Thus we 
can proceed as follows. We begin with pd and go all 
the way to pi imposing that each pj must be uniformly 
distributed in [0,1]. Thus we must have 


9d-l = arcsin i/td-i 


( 8 ) 


and the other angles 9j, with j = I,-- - ,d — 2, should 
be generated as shown in the next equation: 


9j = arcsin ^ 


ti 


UtJj+l C0S2 0fc ' 


(9) 


As was the case with the normalization method, for 
the trigonometric method we need to generate 2{d— 1) 
pRN per pRPV, d — 1 for the angles and d — 1 for 
the random permutation. Nevertheless, because of the 
additional multiplications in Eq. ©, the computation 
time for the last method is in general a little greater 
than that for the former, as is shown in Fig. 01 



where tj, with j = 1, - ■ ■ ,d — 1, are pseudo-random 
numbers with uniform distribution in the interval [0,1]. 
For obvious reasons, this manner of generating a pRPV 
is very unstable, and therefore inappropriate, for nu¬ 
merical implementations. 

A possible way out of this nuisance is simply to ig¬ 
nore the squared cosines in the denominator of Eq. 
That is to say, we may generate the angles using e.g. 

9j = arccos-\/tj (10) 

for all j = 1, • • • ,d — 1. This procedure will give us 
an uniform distribution for cos^ 9j and sin^ 9j , but will 
also increase the chance for the components pj with 


Fig. 4 (color online) Log-log plot of the computation time 
required by the trigonometric method (blue squares) and by 
the normalization method (red circles) to generate a pseudo¬ 
random probability vector with dimension equal to d. The 
calculations were realized using a Processor 1.3 GHz Intel 
Core i5. 


One may wonder if the normalization and trigono¬ 
metric methods, that are at first sight distinct, lead to 
the same probability distributions for the pRPV’s com¬ 
ponents and also if they produce an uniform distribu¬ 
tion for the generated points in the probability hyper¬ 
plane. We provide some evidences for positive answers 
for both questions in Figs. [S] and 01 respectively. 
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Fig. 5 (color online) Semi-log plot of the probability distri¬ 
butions for a component of the unbiased probability vectors 
generated using the trigonometric (lines) and normalization 
(points) methods for some values of d. We see that the two 
methods yield, for all practical purposes, the same probability 
distributions for the components of the pRPV. 



Fig. 6 (color online) Sample with five thousand pRPV gen¬ 
erated using the indicated method. We see that, in this case 
for which d = 3, with exception of the slight overpopulated 
corners, we get a fairly uniform distribution of points in the 
probability space. 

6 Concluding remarks 

In this article we discussed thoroughly three methods 
for generating pseudo-random discrete probability dis¬ 
tributions. We showed that the iid method is not a suit¬ 
able choice for the problem studied here and identified 
some difficulties for the numerical implementation of 
the trigonometric method. The fact that in a direct 


application of both the normalization and trigonomet¬ 
ric methods one shall generate biased probability vec¬ 
tors was emphasized. Then the shuffling of the pseudo¬ 
random probability vector components was shown to 
solve this problem at the cost of the generation of ad¬ 
ditional d — 1 pseudo-random numbers for each pRPV. 

It is worthwhile recalling that pure quantum states 
in can be written in terms of the computational basis 
{|cj)}^=i as follows: 

IV') = T,jCj\cj) = T,j = Ei 

( 11 ) 

where Cj G C and (pj G M. The normalization of |^/>) 
implies that the set {pj} is a probability distribution. 
Thus the results reported in this article are seem to 
have a rather direct application for the generation of 
unbiased pseudo-random state vectors. 

We observe however that the content presented in 
this article can be useful not only for the generation of 
pseudo-random quantum states in quantum informa¬ 
tion science, but also for stochastic numerical simula¬ 
tions in other areas of science. An interesting problem 
for future investigations is with regard to the possibility 
of decreasing the number of pRN, and thus the com¬ 
puter time, required for generating an unbiased pRPV. 
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